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ABSTRACT 

We present a numerical implementation of radiative transfer based on an explicitly 
photon-conserving advection scheme, where radiative fluxes over the cell interfaces of 
a structured or unstructured mesh are calculated with a second-order reconstruction 
of the intensity field. The approach employs a direct discretisation of the radiative 
transfer equation in Boltzmann form with adjustable angular resolution that in prin- 
ciple works equally well in the optically thin and optically thick regimes. In our most 
general formulation of the scheme, the local radiation field is decomposed into a linear 
sum of directional bins of equal solid-angle, tessellating the unit sphere. Each of these 
"cone-fields" is transported independently, with constant intensity as a function of di- 
rection within the cone. Photons propagate at the speed of light (or optionally using a 
reduced speed of light approximation to allow larger timesteps), yielding a fully time- 
dependent solution of the radiative transfer equation that can naturally cope with 
an arbitrary number of sources, as well as with scattering. The method casts sharp 
shadows, subject to the limitations induced by the adopted angular resolution. If the 
number of point sources is small and scattering is unimportant, our implementation 
can alternatively treat each source exactly in angular space, producing shadows whose 
sharpness is only limited by the grid resolution. A third hybrid alternative is to treat 
only a small number of the locally most luminous point sources explicitly, with the 
rest of the radiation intensity followed in a radiative diffusion approximation. We have 
implemented the method in the moving-mesh code AREPO, where it is coupled to the 
hydrodynamics in an operator splitting approach that subcycles the radiative transfer 
alternatingly with the hydrodynamical evolution steps. We also discuss our treatment 
of basic photon sink processes relevant for cosmological reionisation, with a chemical 
network that can accurately deal with non-equilibrium effects. We discuss several tests 
of the new method, including shadowing configurations in two and three dimensions, 
ionised sphere expansion in static and dynamic density field and the ionisation of a 
cosmological density field. The tests agree favourably with analytic expectations and 
results based on other numerical radiative transfer approximations. 
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1 INTRODUCTION 

The transport of radiation and its interaction with matter 
is of fundamental importance in astrophysics, playing a cru- 
cial role in the formation and evolution of objects as di- 
verse as stars, black holes, or galaxies. It would therefore 
be highly desirable to be able to calculate radiative trans- 
fer (RT) processes with equal accuracy and ease as ordinary 
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hydrodynamical and gravitational dynamics. Unfortunately, 
the difficult mathematical structure of the radiative trans- 
fer equation, which takes the form of a partial differential 
equation in six dimensions (3 spatial dimensions, 2 angular 
dimensions, 1 frequency dimension) makes this an extremely 
challenging goal. In fact, the RT problem is so hard, even 
in isolation, that coupled radiation hydrodynamics methods 
are still in their infancy in cosmology thus far. 

However, a large array of different approximations to 
the RT problem have been developed over the years, which 
are often specifically tuned to the requirements and charac- 
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teristics of particular types of problems, and in many cases 
are applied to static density fields only. In this study, we are 
primarily concerned with RT in calculations of cosmological 
reionisation and in star formation, leaving aside other im- 
portant areas such as stellar spectra and atmospheres. Espe- 
cially for the reionisation problem, recent years have seen a 
flurry of activity in the development of new RT solvers that 
are well suited to this problem. These numerical methods 
include long and short characteristics schemes, ray-tracing, 
moment methods and direct solvers, and other particle- or 
Monte-Carlo-based transport methods. In the following we 
briefly describe these schemes. 

In the long characteristics method 
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the computational volume is connected to all other relevant 
cells. Then the RT equation is integrated individually from 
that cell to each of the selected cells. While this method is 
relatively simple and straight-forward, it is also very time 
consuming, since it requires 0{N'^) interactions between 
the cells. Moreover, parallelisation of this approach is 
cumbersome and requires large amounts of data exchange 

between the different processors. 

Short characteristics methods 
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equation of radiative transfer only along lines that connect 
nearby cells, and not to all other cells in the computational 
domain. This reduces the redundancy of the computations 
and makes the scheme easier to parallelise. 

A widely used incarnation of the long-characteristics 
method are so-called ray-tracing schemes. Here a discrete 
number of rays is traced from each source, along which the 
RT equation is integrated in ID, considering absorptions and 
recombinations. As the angular resolution decreases with in- 
creasing distanc e from the source, r ays may be split into 
subrays (e.g. Abel fc Wandeltl [2003 : iTrac fc Cenll2007l ) for 
higher efficiency. The ray-tracing itse l f can be performed 



eithe r on grids (jMellema et al.l I2OO6I : IWhalen fc Norman 



2006 ) or using particles as interpolation points (IBaek et al 



20091 : iGritschneder et alll2009l : [Altav et al.|[20oi ). Other in- 
novative methods trace photons on unstructured grids, for 
example Delaunay tessellations, that are adapted to the 
mean photon optical depth of the gas ([Riikhorst et al.|[200^ : 
iPaardekooper et al.ll2010l : lRitzerveld fc Ickell200^ 

Stochastic integration methods, specifically Monte 
Carlo methods, employ a ray-casting strategy where the 
rays are discretis ed into photon packets |Maselli et al, l2003l : 
iBaek et al.1 12009|) or particles (|Navakshin et al.1 120091) . For 
each photon packet, its frequency and its direction of propa- 
gation are determined by sampling the appropriate distribu- 
tion function of the emitters that have been assigned in the 
initial conditions. A particular advantage of this approach 
is that comparatively few approximations to the radiative 
transfer equations need to be made, so that the quality of 
the results is primarily a function of the number of photon 
packets employed, which can be made larger in proportion 



to the CPU time spent. A disadvantage of these schemes is 
the comparatively high computational cost and the sizable 
level of noise in the simulated radiation field, which only 
slowly diminishes as more photon packets are us ed. The 
'cone' transport scheme of jPawlik fc Schav3l2008l ). where 
radiation is directly transferred between particles, tries to 
improve on these limitations. If needed, this method can 
also create further sampling points dynamically to improve 
the resolution locally. 

Using moments of the radiative transfer equations in- 
stead of the full set of equations can lead to very sub- 
stantial simplifications that can drastically speed up the 
calculations. In this approach, the radiation is represented 
by its mean intensity field throughout the computational 
domain, which is evolved either in a diffusion approxi- 
mation or based on a suitably estimated local Edding- 
ton tensor jGnedin fc Abell I2001I: lAubert fc Tevssieil I2OO8I : 



IPetkova fc Springell l2009l : iFinlator et al.ll2009l ). Instead of 
following rays, the moment equations are solved directly on 
the grid, or in a mesh-less fashion on a set of sampling par- 
ticles. Due to its local nature, the moment approach is com- 
paratively easy to parallelise, but its accuracy is highly prob- 
lem dependent, making it difficult to judge whether the sim- 
plifications employed still provide sufficient accuracy. The 
simplest an d most popular moment method is radiative di f- 
fusion (e.g. IWhitehouse et al.ll2005l : iRevnolds et al.ll2009l ). 
where the RT equation is approximated in terms of an in- 
tegrated energy density in each discretised mass or volume 
element, and this radiation energy density is then evolved 
through the flux-limited diffusion approximation, where the 
flux limiter is introduced to prevent the occurrence of trans- 
fer speeds larger than the speed of light. While the diffusion 
approximation works very well in the optically thick regime, 
its accuracy is hard to judge in general situations. 

A comparison of the relative accuracy of these estab- 
lished RT methods has been carried out i n the "cosmologi- 
cal ra diative transfer comparison project" (|lliev et al.ll2006l . 
I2OO9!), where a subset of the above implementations has been 
compared on a variety of simple test problems. Such a com- 
parison is particularly useful as the lack of known analytic 
solutions for the RT problem, except for a small number 
of simple situations, makes the validation of a RT method 
quite difficult. Reassuringly, the comparison project found in 
most cases reasonably good agreement between the different 
methods, but it also highlighted that each of the different 
techniques shows individual strengths and weaknesses, pro- 
viding ample motivation to search for still better methods. 

It is the goal of this study to propose a new numeri- 
cal scheme for RT that is competitive with the best of the 
known methods in terms of accuracy and general applicabil- 
ity, but is also fast enough to allow self-consistent radiation- 
hydrodynamic simulations in the context of cosmological 
reionisation and star formation problems. We also aim to 
couple the me thod to the new moving-mesh code AREPO 
(|Springelll2010l ) , which solves the equations of hydrodynam- 
ics on an unstructured Voronoi mesh that moves with the 
flow and automatically adapts its resolution to the gravita- 
tional clustering of matter. This mesh-based code computes 
hydrodynamics similar to high-accuracy Eulerian codes on 
Cartesian grids, but it features reduced advection errors 
when the flow velocity is large. 

Our new method is based on a radiation advection tech- 
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nique where a second-order accurate, piece-wise linear re- 
construction of the photon intensity field is used to estimate 
upwind photon fluxes for each face of the mesh. If there is 
only a single point source, such a scheme can exploit the fact 
that the local streaming direction of the photons is known 
everywhere - it is along the ray from the source's position to 
the local coordinate. If there are multiple sources, the radia- 
tion field can be treated equally accurately by decomposing 
it into a linear sum for each source, and treating each com- 
ponent independently. Alternatively, we introduce a direct 
discretisation of angular space, allowing a description of ar- 
bitrary source fields, albeit at the cost of a finite angular res- 
olution. We note that in all these variants the conservation 
of photon number is manifest in the transport step. We treat 
the source terms and the coupling to the hydrodynamics in 
an operator split approach, where the emission, advection, 
and absorption of radiation are calculated in separate steps. 
This makes our approach fully photon conserving, which is 
especially useful for the cosmic reionisation problem, as it 
ensures that all photons emitted by an ionising source are 
really used up in exactly one ionisation event. 

We note that the advection scheme discussed in this 
paper normally propagates the photons at their physical 
speed of light, based on an explicit time integration scheme. 
While this has the advantage of allowing general, fully 
time-dependent radiative transfer simulations, it can make 
them computationally very expensive due to the required 
small Courant time steps. This can however be greatly al- 
leviatcd by using a reduced speed of light approximation 
ijGnedin fc Abelll200li l. which allows much larger timesteps 
while still preserving the speed of cosmological ionisation 
fronts (I- fronts). With this approximation, it then becomes 
possible to calculate high-resolution cosmological radiation 
hydrodynamics simulations of structure formation that si- 
multaneously account for cosmic reionisation, with no re- 
striction on the number of sources. 

In Section [2] of this paper, we present our methodol- 
ogy in detail. We first give a brief introduction to the RT 
equation in Section 12.11 Then we discuss three variants of 
our solution method for the radiation advection equation in 
Sections 12. 2112.31 and 12.41 In Section [2.51 we briefly describe 
our treatment of emission and absorption processes, with 
an emphasis on the hydrogen chemistry relevant for the cos- 
mic reionisation problem, and in Section [2.61 we specify our 
formulation of photo-heating and radiative cooling. Issues 
of time stepping and code implementation are discussed in 
Sections 12.71 and 12.81 We move on to a presentation of basic 
test results in Section[3l starting with a variety of shadowing 
(Section 13. ip and Stromgren sphere tests f Section 13. 2|) . We 
then consider the more demanding tests of I-front trapping 
in Section 13.31 the ionisation of a cosmological density field 
in Section [331 and an ionisation problem with dynamic den- 
sity field in Section [3.51 Finally, we present our conclusions 
in Section ID 



2 AN ADVECTION SOLVER FOR THE 
RADIATIVE TRANSFER PROBLEM 

2.1 The radiative transfer equation 

Let us briefly discuss different forms of the RT equation, 
which is helpful to clarify how our new method differs from 



other approaches, and for specifying our notation. Let f~f = 
ff{t, X, p) be the photon distribution function for comoving 
coordinate x and photon momentum 



p = a — n , 

c 



(1) 



where a = a{t) is the cosmological scale factor, h is the 
Planck constant, v is the frequency of the photons, and n is 
the unit vector in the direction of photon propagation. Then 
the number of photons in some part of the Universe is 



A'^ — I dx dp /t- {t, X, p) 



(2) 



We can quite generally write the phase-space continuity 
equation for the distribution function of photons as 



dt 



+ 



d_ 
9x 



x/.) + ^(p,f.)=^ 



dt 



(3) 



In this Boltzmann-like transport equation, the source and 
sink terms on the right-hand side of the equation represent 
photon emission and absorption processes, respectively. If 
we neglect gravitational lensing effects, individual photons 
propagate along straight lines with conserved momenta, i.e. 
we have x = (c/a)n and p = 0. The transport equation 
hence simplifies to 



dt 



c 9 ^ , df. 
+ - — (n^) = 



a 9x 



dt 



sources 



dt 



(4) 



Normally, a direct use of equation @ through a dis- 
cretisation of phase-space is considered prohibitively expen- 
sive due to the high-dimensionality of the problem. How- 
ever, if only monochromatic radiation is considered, which 
is often sufficient, the momentum-space dimensions reduce 
to just two angular coordinates. If furthermore only a rela- 
tively coarse angular resolution for the photon transport is 
sufficient, then the 4n solid angle described by these angular 
dimensions may be discretised into a limited set of cones, say 
up to 10-100, at which point a brute-force solution of equa- 
tion Q on a 3D mesh becomes computationally feasible and 
attractive, as we shall argue here. 

Before we discuss this in more detail, let us first briefly 
recall for clarity how the specific intensity that is normally 
used in RT studies relates to equation Q. We can define 
the specific radiation intensity in a certain direction fi 
through the energy AE^ = I^AvAAAilAt of photons that 
pass through a physical area normal to n and within 
solid angle AQ around fi, over a time interval At and in a 
frequency bin Ai/. With this definition, the specific intensity 
Ii, is then related to the photon distribution function as 



Ii, — hufj 



f-r 



(5) 



du dn dA dt 

Substituting into equation @ , and writing the absorp- 
tion and emission terms in their conventional form, one ob- 
tains the cosmological RT equation in the form 



1 dli, n dli, 
c dt a 9x 



H{a) ( dh \ 



(6) 



where is the absorption coefficient, is the emission 
coefficient, and H{a) is the Hubble rate. Defining the solid 
angle averaged intensity as 



1 

47r 



dnh, 



(7) 
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Cell 2 




Figure 1. A simple sketcli sliowing tlie geometry involved in our 
advection scheme for a single point source at coordinate Xs. Here 
n is the photon propagation direction, f is the normal vector of 
a face of a cell, Xc is the center of mass of the corresponding cell 
and is the center of mass of the face for which the photon flux 
is calculated. 



we can calculate the physical photon number density from 
the specific intensity as 



-.phys 



hv 



Av. 



(8) 



Another equivalent way to obtain the photon number den- 
sity is simply to integrate the distribution function, 



dp/^(f,x,p) 



(9) 



which yields the comoving number density of photons, = 
a^'nJi^^'' . This highlights again that describing the radiation 
field with the arguably more familiar RT equation ((6)1, or 
with the distribution function and the Boltzmann-like equa- 
tion is fully equivalent. In this paper, we will mostly 
work in the latter formulation. 

In general, to solve the RT problem on some discretised 
mesh, we can split off the source and sink terms and treat 
them separately in the time integration. In such an operator 
splitting approach, known as Strang splitting, we are basi- 
cally left with two separate problems that are interleaved in 
the time integration, one is to follow the conservative trans- 
port of photons on the mesh, the other is the local updat- 
ing of the photon density field through the source and sink 
terms. In the following, we first focus on the conservative 
transport problem, which is where the primary computa- 
tional challenge lies. 

2.2 Transferring radiation by advection for point 

sources 

Suppose for the moment that we know at a given point in 
space that all photons stream in the same direction fi. This 



is for example the case if there is a single point source at 
coordinate Xs (i.e. no other sources and no scattering are 
present). For simplicity, we shall also restrict ourselves to a 
spatially invariant photon momentum spectrum. One then 
obtains a simple advection equation for the comoving photon 
density n-y-. 



dn. 



cn „ 

H Vn^ 

at a 



0, 



(10) 



where the local advection direction n is known at every point 
X and is simply given by 



n(x) = 



(11) 



This advection is conservative and may be solved with the 
techniques commonly employed to treat the hyperbolic con- 
servation laws of ideal fluid dynamics on spatial meshes. In- 
deed, this is the approach we are going to employ: we shall 
use a conservative transport scheme based on a second-order 
accurate upwind method that is inlined with the hydrody- 
namic calculations of our unstructured moving-mesh hydro- 
dynamics code AREPO, which is described in some more 
detail below. It is important to note that knowledge of the 
local number density field of photons combined with the 
source location Xa is sufficient to accurately solve the ra- 
diative transport, simply because this information suffices 
to specify the photon streaming direction at every point in 
space. Apart from the spatial discretisation, no approxima- 
tions need to be made for the case of a single monochromatic 
point source in this treatment. 

In practice, we use a second-order accurate spatial re- 
construction technique to convert photon numbers A^^ stored 
for each cell i of a given mesh into a photon density field. 
For every cell, we first obtain an estimate {Vn-,)^ for the 
gradient of n^, which allows a piece- wise linear conservative 
reconstruction of the photon density field, in the form 



n-f{x) = (n^.) . -I- {Vn^)^ (x — x^), for x £ cell i. 



(12) 



Here {n^)^ = Ni/Vi is the mean photon number density 
in the cell with center-of-mass x^ and volume Vi. As il- 
lustrated in the sketch of Figure [T] for every face centroid 
X/ of the mesh, we can then identify the upwind side of 
the photon flow, based on the sign of the dot product be- 
tween face normal f and the photon streaming direction 
n — (x/ — Xs)/|x/ — Xs|. This allows us to estimate the 
photon flux F-y over the face as 



7 - - (f -n) n^i^f) 
a 



(13) 



where the photon density ^^.(x/) at the face centroid is es- 
timated based on the linear reconstruction of the cell on the 
upwind side. If the face has comoving area A, the number 
of photons exchanged during time At between the cells that 
share the face is then given by 



AN^ = F^AAt. 



(14) 



Due to the pairwise exchange of photons, the conservation 
of total photon number is manifest, which is important 
for guaranteeing that I-fronts propagate at their physical 
speeds. We note that in our code the mesh is composed of 
Voronoi cells (of which a Cartesian mesh is a special case), 
but this is not important for the general approach. 
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^ource 1 



Figure 2. Sketch that illustrates the linear summation principle 
used to treat the radiative transfer for multiple sources. Here 
ni and ri2 are the photon propagation directions from the two 
sources, as seen from the center of a face with normal vector 
f . The total flux passing through the face is then computed as a 
linear sum of the contributions from the partial fields created by 
each source. 



There are two important caveats with this transport 
scheme, which need to be pointed out. One is that this ex- 
phcit transport scheme requires a time step that is given by 
a local Courant criterion for the photons, which can become 
very small due to the high speed of light. For reionisation 
problems, this can however be circumvented with the re- 
duced speed of Ught approximation, which we will discuss 
in more detail later on. The other caveat is that close to a 
point source the mesh resolution will always be coarse, so 
that our use of a single Gauss point per mesh face may intro- 
duce sizable errors in the discretised advection fluxes. This 
can happen when the opening angle under which a mesh 
face is seen by the point source is large, so that adopting a 
single propagation direction for the entire face is inaccurate. 
As a result, isophots of the radiation field produced by the 
point source may then deviate from sphericity with distor- 
tions that reflect the local geometry of the mesh around the 
point source. We have found however that this problem can 
be cured quite effectively by injecting the photons of the 
source in a kernel-weighted fashion over 2-3 mesh cells or 
so. With such slightly extended sources, the above scheme 
is able to quite accurately treat single point sources. 

The approach can also be straightforwardly extended to 
multiple point sources simply by linear superposition of the 
radiation fields produced by each of the individual sources, 
as sketched in Figure [2] This means that the advection equa- 
tion is solved for the radiation field of each source separately. 
This obviously involves a computational and storage cost 
that scales with the number of sources, but if the number of 
sources is small, this is an interesting technique for certain 
applications due to its high accuracy. As we show in our test 
problems, the method in particular is able to accurately cast 



shadows, and unlike for example in the optically thin vari- 
able Eddington tensor approximation (OTVET) , there is no 
accuracy-degrading mutual influence of multiple sources on 
each other. 

However, for a large number of ionising sources, the lin- 
ear superposition approach will quickly become infeasible. 
For example, in large cosmological simulations, we would 
like to allow every star particle to act as a source of ionising 
radiation. Here we obviously cannot decompose the radia- 
tion field into all its single point sources, instead, we need 
to employ another decomposition. We have actually devel- 
oped two possible schemes for this, which we describe in the 
following. 



2.3 A hybrid between point-source treatment and 
local diffusion 

One possibility to address the multiple point sources prob- 
lem is to only retain a finite number A'br of locally bright- 
est sources in an explicit treatment, while all the remaining 
sources are lumped together into a background radiation 
field that is treated with radiative diffusion. The idea here 
is that especially in cosmic reionisation problems the local 
ionisation "bubble" is expected to be driven primarily by 
one or a few sources, and only at very late stages, multiple 
sources may become visible at a given point, but then reioni- 
sation has largely completed already anyway. By making the 
set of sources that are treated exactly as point sources spa- 
tially variable, we should then get a quite accurate approx- 
imation of the reionisation phenomenon even for moderate 
values of A'^br • Since in the limit of large Abr , the scheme will 
become essentially exact (apart from spatial discretisation 
errors), the degree to which imposing a limiting value for 
Abr affects the results can be readily tested. 

We will discuss some results obtained with this ap- 
proach later on, but we note that it clearly involves sev- 
eral complication when applied in practice. First of all, the 
need to allow a local change of the list of bright sources re- 
quires that one keeps track of all locally incoming fluxes of 
radiation, sorting them appropriately, and matching them 
through the use of unique source identifiers to the already 
stored radiation intensities from the previous step. Also, 
since neighboring cells may have different source lists, a 
matching procedure is required for gradient estimates, with 
the additional complication that the accuracy of the gradi- 
ents will be reduced at "domain boundaries", i.e. regions 
of the mesh that differ in their assessment what the lo- 
cally most important Abr point sources are. Furthermore, 
if the number of sources is very large and spread out in 
space (e.g. the individual stars in galaxies), the injection of 
photons needs to be treated in some sort of clustered fash- 
ion, otherwise faint individual sources may not be able to 
compete with the Abr bright sources already stored locally, 
so that they are channeled into the radiatively treated fiux 
reservoir right away without having a chance to build up to a 
significant source when combined with the potentially many 
nearby sources that are equally faint. Finally, one also needs 
a separate radiative diffusion solver, which requires a small 
timestep for stability when integrated explicitly in time, as 
we do here. For all these reasons, we actually favour in most 
applications our second approach for treating a large num- 
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ber of sources, which is faciUtated by discretising the solid 
angle explicitly, as we describe next. 

2.4 Full angular discretisation and cone transport 

For general radiation fields we seek a method that can di- 
rectly represent the angular distribution of the local radia- 
tion field. This can, for example, be done in terms of mo- 
ments of the radiation field. However, we here want to pro- 
pose a more flexible approach that is based on a direct an- 
gular discretisation of the photon space. To this extent, we 
can decompose the full solid angle into a set of cones of 
equal size, for exam ple based on the well-known HEALPIX 
l|G6rski et al.l |2005| ) tessellation of the unit-sphere, which 
we shall use in the following. Our strategy could however be 
straightforwardly generalised also for other discretisations of 
angular space. In HEALPIX, the unit sphere is decomposed 
into A'^pix — 12 A'^de patches of equal solid angle (which we 
call "cones" for simplicity, even though they are not exactly 
axi-symmetric), each centred around a central direction Aj, 
where j = 1 . . . A^pix. We now linearly decompose the ra- 
diation field /-y into A^pix components, each containing the 
photons that propagate along a direction within the corre- 
sponding cone: 

^(x,n) = ^/;^(x,n), (15) 

j 

where /;^(x,n) = if the photon direction n lies outside 
of Ai}j around n^. The basic simplification we now make 
is that we assume that each of the partial radiation fields, 
/;^(x, n), can be taken to be constant as a function of di- 
rection within the corresponding cone. Or in other words, 
each of the partial fields (x, n) describes the intensity of 
a homogeneously illuminated beam of opening angle AQj 
around direction fij , emanating from the local coordinate x. 
Our goal is now to generalise the radiation advection scheme 
for point sources outlined above such that it can accurately 
transport the radiation cones occurring in this discretisa- 
tion. 

If we simply transport one of the partial radiation fields 
f^i locally always along the primary direction of its cone, i.e. 

f + ^.V/-0, (16) 

we will invariably observe a central "focusing effect" , i.e. the 
radiation emanating from a point will not illuminate the fi- 
nite solid angle AQj homogeneously, but rather tend to con- 
centrate along the primary axis of the cone. It is clear that 
this "focusing effect" arises from the parallel transport de- 
scribed by equation (|16|l : instead of transporting the photon 
field over different directions that are uniformly spread over 
the finite solid angle, all of the photons are transported along 
the single direction nj, with any residual angular spread 
around nj arising only from numerical diffusion due to the 
finite mesh resolution. 

One may try to fix this problem by somehow randomis- 
ing the direction within the corresponding cone taken in 
single transport steps, or by using higher-order quadratures 
in integrating the fiuxes arising for a given mesh geometry. 
However, we have found that a simple trick can be used to 
resolve this issue, and to obtain close to perfect results even 




Figure 3. This sketch illustrates the geometry and the vectors 
involved in our "cone transport" , where the angular space is dis- 
cretised into regions of equal angle (in 2D) or solid angle (in 3D). 
In this example, only four cones in in 2D are used. The photon 
field is linearly decomposed into radiation fields corresponding 
to the four cones, which have symmetry axes fii, n2, na, ... iijv, 
where N is the number of discrete cones or angles, i.e. = 4 
in the sketch. At each face of the mesh (here the normal vectors 
fi and f2 are shown), photon fluxes for each of the partial fields 
are estimated. The photon propagation direction is taken to be 
parallel to the gradient of the total radiation intensity field, con- 
strained to lie within the opening angle of the corresponding cone. 



n' 




Figure 4. A sketch illustrating the construction of the vector 
given by equation II20I I. The symmetry axis of the solid-angle cone 
j is given by Aj , while the gradient direction is n'j . If the latter lies 
outside the cone, it is projected onto the cone to yield direction 
(iijOnew, which is then used in the local advection step for the 
cone's radiation field. 
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for unfavourable mesh geometries. To this end, we replace 
the local advection direction lij appearing in equation (|16|) 
with a modified direction n^ , chosen along the gradient of 
the tota/ radiation density field, but constrained to lie within 
the cone j. Specifically, we first adopt 



V/7 



(17) 



and calculate the angle between the gradient direction and 
the cone direction as 

4> — arccos (n^ • fij). (18) 

If this angle is larger than the half opening angle of the cone, 

= v/(47r/iVpix)A, (19) 

then we use the vector (nj )new, which is defined by the inter- 
section of the plane spanned by n^- and fij with the cone of 
half-opening angle (jf"-'^^ (see Figure!?}. This vector is given 
by 



(nj)ncw = sin {4>^'^'^)va. + cos ( 
where 

g = fij X fij , 
m = g X fij. 



(20) 

(21) 
(22) 



In other words, we transport the radiation corresponding 
to a certain cone always in the direction of the negative 
intensity gradient, constrained to lie within the solid angle 
defined by the cone. 

It is clear that this modification has the tendency to 
smooth out the angular gradient of the radiation field within 
a cone, making it uniform in the cone. For example, imag- 
ine that the transport has led to some intensity excess along 
the principal direction of the cone. This will then cause some 
of the transport steps to propagate photons away from the 
symmetry axis of the cone, slightly more sideways, until the 
cone is illuminated homogeneously again. But importantly, 
the constraint we imposed on the advection direction means 
that all of the photons of any of the partial radiation fields 
are always transported along a direction "permitted" by 
their corresponding angular cone. While the specific choice 
for this direction may hence deviate slightly from the pri- 
mary cone axis n^, this deviation is strictly bounded, and 
it will automatically become smaller if a larger number of 
angular cones is used. One may wonder why we base the ini- 
tial calculation of the transport direction in equation \n\ 
on the total radiation intensity field, and not on the partial 
cone field fi/ alone. This is done to avoid possible boundary 
effects at the edges of cones, for example when two neigh- 
bouring cones are both homogeneously illuminated. Using 
the gradient of the total field will in this case automatically 
work to eliminate any residuals from the common boundary 
and to produce a seamless connection of the cones, a fea- 
ture that is not guaranteed when the gradient of the partial 
field is used instead. In Section (3] we will discuss a number 
of test problems that illustrate that our simple approach 
works rather well in practice. 

We note that the angular discretisation we outlined here 
is completely independent of the total number of sources. 
Also, its angular resolution is constant everywhere (at least 
in the present implementation), even though the spatial res- 



olution of the mesh can vary as a function of position. An- 
other interesting aspect of the method is that it can work 
accurately both in the optically thin and in the optically 
thick regime, as well as in the transition region. Unlike in 
certain approximate treatments of RT, for example in ra- 
diative diffusion, we have not made any approximation that 
changes the fundamental character of the equations, apart 
from the use of a spatial and an angular discretisation. This 
suggests that the robustness and the convergence of results 
obtained with this method can reliably be tested by simply 
changing the grid and/or angular resolution, and if conver- 
gence is achieved, then the method should converge to the 
correct solution in the limit of high resolution. The latter 
property can not necessarily be expected for RT schemes 
that use more drastic approximations. 



2.5 Source and sink terms 

We treat source and sink terms in the radiative transfer 
equation through an operator splitting approach, where the 
evolution of the homogeneous RT equation (which conserves 
photon number) is alternated with an evolution of the source 
terms alone. This greatly simplifies the calculation of the 
interaction of the local radiation field with matter, and also 
allows accurate balance equations that for example ensure 
that the number of photons absorbed matches the number of 
atoms that are ionised. As an illustrative example, we here 
detail our implementation of hydrogen chemistry, which can 
be used in simple model calculations of cosmic reionisation. 



2.5.1 Emission processes 

Emission of ionising radiation in a cosmological simulation 
can be based on a variety of source models, tied for ex- 
ample to star-forming gaseous cells, star particles, or sink 
particles that represent accreting supermassive black holes. 
Given the source luminosities and their coordinates, we can 
simply find the cells in which the sources fall, and inject the 
number of photons emitted by them over the timestep into 
each of the corresponding host cells. We normally assume 
isotropic sources where we distribute the total emissivity 
equally over all angular cones, but in principle also beamed 
emission characteristics can be realised. 

If our single/multiple point source approach is used in- 
stead, we spread the source photons over a small region 
around the host cell with a Gaussian-shaped kernel with 
a radius equal to a few effective host-cell radii. This is done 
to avoid potential asymmetries in the source's radiation field 
that otherwise can arise from the particular geometry of the 
source cell. 



2.5.2 Absorption and hydrogen chemistry 

For simplicity, we here discuss a minimal chemical model 
that only follows hydrogen and an ionising photon density 
field with a fixed spectral shape. Extensions to include he- 
lium and several ionising frequencies to account for changes 
of the spectral shape can be constructed in similar ways. 

The neutral hydrogen fraction nni evolves due to photo- 
ionisations, coUisional ionisations and recombinations: 
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dhm 
dt 



— anu ficnHii ~ finu hchm — canu fimhj, 



(23) 



where a{T) is the recombination coefficient, /3(T) is the col- 
lisional ionisation coefficient and a is the effective photo- 
ionisation cross-section of neutral hydrogen for our adopted 
spectrum, defined as 



hv 



)Av 



hv 



■Av 



(24) 



Here a^iy) is the frequency dependent photo- ionisation 
cross-section of neutral hydrogen (with cr^ = for frequen- 
cies u < below the ionisation cut-off vo). The photon 
density on the other hand evolves according to 

dn-y _ _ 

= -canunum^- (25) 

Ifere the variables nni, nnii, and h^ express the cor- 
responding abundance quantities in dimensionless form, in 
units of the total hydrogen number density tih, for example 
/riH- If we consider only hydrogen, we hence have 
the constraints hg = whii and hui + mhii = 1- 

In order to robustly, efficiently and accurately integrate 
these stiff differential equations, special care must be taken. 
This is especially important if one wants to obtain the cor- 
rect post-ionisation temperatures, which requires an accu- 
rate treatment of the rapid non-equilibrium effects during 
the transition from the neutral to the ionised state (e.g. 
iBolton et al]|2005h . Also, one would like to ensure that the 
number of photons consumed matches the number of hydro- 
gen photo-ionisations, and that the injected photo-heating 
energy is strictly proportional to the number of photons ab- 
sorbed. We use either an explicit, semi-implicit, or exact in- 
tegration of equations 1)23^ and (|25|) to achieve these goals, 
depending on the current conditions encountered in each 
step. 

Specifically, we start by first calculating an explicit 
estimate of the photon abundance change over the next 
timestep, as 



hL — —canu huihL At, 



(26) 

where i enumerates the individual timesteps. If the implied 
relative photon density change is small, say \ Ah^ \ < 0.05 hi,, 
we are either in approximate photo-ionisation equilibrium or 
the photon density is so large that it does not change appre- 
ciably due to hydrogen ionisation losses during the step. In 
this situation, we can calculate an estimate for the neutral 
hydrogen density at the end of the step based on implicitly 
solving 



nui + [a(l-n 



i + l\2 
HI ) ~ 



-l3h]+\l-h'+j^)]nHAt+Ah^{27) 



for n^^. If the implied relative change in hlu is again small, 
we keep the solution. 

Otherwise, we first check whether the photon number 
is very much smaller than the neutral hydrogen number, 
i.e. whether we have n-y < 0.01 tt-hi. If this holds, the photons 
in the cell cannot possibly ionise a significant fraction of the 
neutral hydrogen atoms, but the photon abundance itself 
may still change strongly over the step (for example because 
almost all of the photons are absorbed). We in this case first 
compute an estimate of the new photon number at the end 
of the step, based on the implicit step 



With the solution for hl^'^ in hand, we calculate again an 
implicit solution for the new neutral hydrogen fraction at 
the end of the step, using equation (|27p . If the predicted 
relative change in the hydrogen ionisation state is small, we 
keep the solution, otherwise we discard it. 

Finally, if both of the two approaches to calculate new 
values for hl^^ and h'^^ at the end of the step have failed, 
we integrate the rate equations ((23} and (|25p essentially ex- 
actly over the timestep At, using a 4-th order Runge-Kutta- 
Fehlberg integrator with adaptive step-size control as imple- 
mented in the GSL librar}0. We note that this sub-cycled 
integration is hence only done in timesteps where the ion- 
isation state changes rapidly in time and non-equilibrium 
effects can become important, which is a very small fraction 
of all cells, such that our updating scheme remains compu- 
tationally very efficient. 



2.6 Photo-heating and radiative cooling 

To calculate the evolution of the thermal energy, we can now 
inject the photo-heating energy 



AEj = {hi, 



)nHV 



(29) 



into the corresponding cell, where V is the volume of the cell 
under consideration, {hi, — hl^^)n-ii is the number density of 
photons consumed by ionising events over the timestep, and 
e-y gives the average energy absorbed per photo-ionisation 
event. For our prescribed spectral shape, this injection en- 
ergy per ionisation eve nt is given by the frequency-averaged 
photon excess energy (|Spitzerlll998l ) 



, 47rJ^ 

du— (Ti 

nv 



{hv — huo) 



, 4nJ^ 
dv— 

nv 



(30) 



above the ionisation cut-off vq. For many of our test calcula- 
tions, we assume a black body spectrum with T^a — 10^ K, 



which leads to e~, 



6.4 eV. 



hi, — cannhlijnl^^ 



At. 



(28) 



The evolution of the thermal energy is then completed 
by a separate cooling step that accounts for recombina- 
tion cooling, coUisional ionis ation, excitation cooling, and 
bremsstrahlung cooling (e.g. iKatz et al.lll99d ). We imple- 
ment these cooling rates with a combination of an explicit 
and implicit timestep integrator, where an explicit integra- 
tion scheme is used as default, but if the temperature change 
over the step becomes large, the cooling is instead calculated 
with an unconditionally stable implicit solver. 



2.7 Time stepping and the reduced speed-of-light 
approximation 

As discussed above, we include the source terms into the 
time integration of our RT solver by an operator splitting 
technique, where the source and advection parts are treated 
separately. This technique can be generalised also to the 
coupling of hydrodynamics and radiative transfer, by alter- 
natingly evolving the hydrodynamical density field and the 
radiation field with its associated radiation chemistry. In 
fact, this is the approach we follow in our radiative trans- 
fer implementation in the hydrodynamical AREPO code. As 



http://www.gnu.org/software/gsI 
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the latter is a moving-mesh code, we however need to ensure 
that during the hydrodynamical step the radiation field is 
left invariant. This can be achieved by appropriate advec- 
tion terms that compensate for the mesh-motion during the 
hydrodynamical step. 

For the time integration of the radiative source terms, 
we employ implicit or semi-implicit methods, as described 
in Section [231 that are stable even for very large time steps, 
and in selected situations, adaptive numerical integration 
of the stiff ordinary differential equations that describe the 
chemical networks. The latter is essential to accurately ac- 
count for non-equilibrium effects. As these processes are 
completely local, this does usually not incur a very signif- 
icant computational cost, provided the exact integration is 
only done where really needed. In contrast, the timestepping 
of the radiation advection step poses more severe computa- 
tional requirements. This is because this step is based on an 
explicit time integration scheme, whose timestep needs to 
obey a Courant criterion of the form 

Afadvoct < Ck , (31) 

c 

where < Cfe < 1 is the Courant factor, and Aa; is taken to 
the smallest comoving size of a cell in the simulated volume. 

In ordinary hydrodynamics, a similar time step con- 
straint is encountered, except that the speed of light is re- 
placed with the speed of sound. Since we are primarily inter- 
ested in non-relativistic gas dynamics in cosmological struc- 
ture formation, the speed of light will typically be a factor 
~ 10^ to 10'' larger than the hydrodynamical sound speed. 
The resulting reduction in the allowed time step size can 
hence make a simulation prohibitively expensive when the 
RT is coupled to the hydrodynamics over significant frac- 
tions of the Hubble time. However, in many applications 
of interest this problem can be greatly alleviated by re- 
sorting to an artificially reduced speed of light c', which 
is introduced instead of the physical speed of light both 
in the transpor t equa tion a nd the ionisation e q uatio n. As 
iGnedin fc Abel (|200ll ) and lAubert fc TevssieJ (|2008l ') dis- 
cuss in detail, this reduced speed-of-light approximation is 
especially attractive for cosmic reionisation problems be- 
cause it here does not modify the propagation speed of I- 
fronts, except perhaps in the very near field region around a 
source directly after it turns on, but this introduces a neg- 
ligible timing error. In general, the reduced speed-of-light 
approximation can be expected to yield reasonable accu- 
racy in many radiation hydrodynamic problems as long as 
c' remains significant larger than the maximum sound speed 
occurring in the simulation. 

2.8 Implementation aspects in the moving- mesh 
code AREPO 

We have implemented the different variants of our radi- 
ation advection solver in the moving-mesh code AREPO 
ISpringejlioiol '). This code treats hydrodynamics with an 
ordinary finite-volume approach and a second order accu- 
rate Godunov scheme, similar to many Eulerian grid codes. 
However, AREPO works on an unstructured mesh created 
with a tessellation technique. The particular mesh used is 
the Voronoi tessellation created by a set of mesh-generating 
points. Using such a mesh offers a number of advantages 



compared to traditional grid codes in that its mesh can flow 
along with the gas. As a result of the induced dynamic mesh 
motion, AREPO exhibits considerably lower advection errors 
than ordinary mesh codes, and also avoids the introduction 
of preferred spatial directions. Also, the cell size automati- 
cally and continuously adjusts to the density in a Lagrangian 
sense, and is hence decreased in regions where typically more 
resolution is required even without doing adaptive mesh re- 
flnement. 

For implementing our RT transfer scheme as described 
above, we can readily employ the infrastructure and com- 
munication algorithms provided by the fully parallelised 
AREPO code, making it an ideal base for a first demonstra- 
tion of the method. This in particular applies to the gradient 
estimation, the spatial reconstruction of the photon intensity 
fields, and the parallelisation for distributed memory com- 
puters. A full des cription of thes e aspects of our code can 
hence be found in ISpringell (|2010l ). We carry out a RT step 
on every top-level synchronisation point of the AREPO code, 
which means on the longest time step Afmax allowed by the 
gravitational and hydrodynamical interactions followed by 
the code. If Afadvoct is smaller than the top-level simulation 
time-step Afmax, the radiation transfer step is calculated in 
several sub-cycling steps equal to or smaller than Afadvcct, 
as needed. 

Note that these sub-cycling steps do not require a new 
construction of the Voronoi mesh, or a new gravity calcu- 
lation, hence they are in principle quite fast compared to 
a full step of the hydrodynamic code. However, this advan- 
tage can be quickly (over)compensated by the need to carry 
out multiple fiux calculations for each of the angular com- 
ponents of the radiation field, and the additional need to do 
subcycling in time to ensure stability of the explicit time in- 
tegration used in the advection steps. Furthermore, if a mul- 
tiple frequency treatment is desired, the cost of the radia- 
tive transfer calculations will scale linearly with the number 
of frequency bins employed, simply because the dominating 
advection part of the radiative transfer problem needs to 
be carried out for each frequency independently. The addi- 
tional storage requirements for a multiple frequency treat- 
ment should also not be overlooked, which again scale lin- 
early with the number of frequency bins, likewise with the 
number of solid-angle bins used in the angular discretisation. 
It is hence clear that multi-frequency radiative transfer at 
high angular resolution clearly remains expensive with the 
discretisation scheme proposed here. However, the relative 
cost increase compared to hydrodynamics alone is a constant 
(and at least for sufficiently interesting problems still afford- 
able) factor that is nearly independent of spatial resolution. 
This, together with the ability of our scheme to cope with es- 
sentially arbitrary source functions, makes it an interesting 
new technique for cosmological hydrodynamics. 



3 BASIC TEST PROBLEMS 

3.1 Shadows around isolated and multiple point 
sources 

We begin our investigation of the accuracy of our proposed 
radiative transfer algorithms with isolated point sources in 
an optically thin medium that includes some regions with 
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Figure 5. Photon density maps of 2D shadowing tests in three different cases; single point source with a single obstacle (left panel), 
two point sources with two obstacles (middle panel), and a single source with two obstacles (right panel). The green lines indicate the 
geometric boundaries of the expected shadow regions, whereas the thick white lines mark the absorbing obstacles. 



absorbing obstacles. This serves both as a verification that 
an isolated point source produces a radiation field cc 1 /r^ 
(in 3D, and oc l/r in 2D), with sufficiently spherical 
isophots, and as a test whether the method can cast sharp 
shadows behind obstacles. The latter is often difficult for 
RT transfer schemes, especially the ones that are diffusiv e in 
character such as the OTVET scheme (e.g. iGn cdin & Abel 
l200ll :[ Petkova fc SpringelllioogI ). 

In Figure [5l we show such shadowing tests for three dif- 
ferent cases, which for visualisation purposes have been done 
in 2D space. In the left hand panel, we consider the shadow 
that is produced by an obstacle when it is illuminated by 
a single source in the middle of the panel. The green lines 
show the geometric boundaries of the theoretically expected 
position of the shadow. We see that the obstacle produces 
a rather sharply defined shadow with only a small radiation 
leak into the shadowed region due to numerical advection 
and discretisation errors along the shadow boundaries. In 
the unshadowed regions, the radiation intensity falls of as 
oc l/r, as expected. 

Equally good results are also obtained when multiple 
sources are considered in our "linear sum" approach to the 
total radiation field, where the total photon density is com- 
puted as a linear sum of the photon fields from each source, 
and the transport of each partial field is treated indepen- 
dently. Examples for this are shown in the middle and right 
panels of Figure [SI where two obstacles and one or two 
sources are used in different configurations. Again, the shad- 
ows agree very well with the expected boundaries shown 
with green lines, with only a small amount of residual diffu- 
sion into the shadowed regions. If the spatial mesh resolution 
is improved, the shadows become progressively sharper still. 

We note that the above success essentially holds in this 
approach for an arbitrary set of absorbing regions, and an 
arbitrary combination of point sources. It hence provides a 
general and highly accurate solution to the radiative trans- 
fer problem, even though it can certainly get expensive to 
obtain it, especially for a large number of source. It is impor- 
tant to note however that the radiation fields produced by 
our scheme are essentially noise- free, which is a drastic im- 
provement compared to results obtained from schemes that 
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Figure 6. Radiation field around two point sources and two ab- 
sorbing obstacles, for our hybrid treatment of point-sources and 
radiative diffusion. In this example, only the brightest source seen 
from a given cell was treated explicitly as a point source, while the 
other radiation was dumped into a background field transported 
with radiative diffusion. 



rely on Monte-Carlo methods (e.g. Maselli et ahl 20031 ) 



on randomised cone transport (|Pawlik fc Schavell2008l ') 

As we discussed earlier, for many problems in astro- 
physics the number of sources is too large to make the lin- 
ear sum approach a viable solution technique. In our first 
approach to work around this limitation, we only treat the 
photons from the brightest sources at a given cells as in- 
dependent point sources in the transport scheme, while all 
other incoming photons form fainter sources are added to a 
background radiation field, which is then diffused from cell 
to cell. In Figure (6] we show a (somewhat extreme) exam- 
ple for how this can change the results. We repeat the test 
shown in the middle panel of Fig. [5] which has two sources 
and two obstacles, but this time we only allow the code to 
treat the locally brightest A^br = 1 sources as explicit point 
sources, while the rest of the radiation needs to treated with 
radiative diffusion. As we can see from Figure (6] the radia- 
tion field near to the two sources is unchanged, as expected, 
but at the mid-plane, where the sources have equal intensity. 
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Figure 7. Illustration of our cone transport scheme, and the ac- 
curacy with which it can represent homogeneously illuminated 
radiation cones. The panel shows a map of the 2D photon den- 
sity around a single source positioned at the center of the field. 
Radiation transport was calculated by partitioning the full 2-k 
angle into eight regions of size 7r/4 each, with each of the fields 
transported individually. To show that the cones produce a homo- 
geneous field, the source luminosity was only injected into every 
second cone and some of the cone boundaries were marked with 
green lines. 



half of the flux is dumped into a diffusive reservoir. The dif- 
fusion approximation then lets the radiation spread from the 
mid-plane more slowly, causing an incorrect increase of the 
radiation intensity there. A second effect is that the shadows 
behind the obstacles are not sustained as nicely any more, 
instead they are partially illuminated by the radiative dif- 
fusion. It is important to be aware that unlike in the pure 
transport scheme considered earlier, these errors will not be- 
come smaller for an improved mesh resolution, rather, one 
would simply converge to a wrong solution in this case. 

The example studied in Fig. [6] is deliberately extreme 
in the sense that A^br was kept very low. Much better re- 
sults can be expected for a sizable value of Nhi, say 5 — 10, 
because then the flux that needs to be treated with the dif- 
fusion approximation should become locally sub-dominant 
everywhere. Nevertheless, for general radiation fields and 
smoothly distributed source functions, we prefer our "cone 
transport" scheme, which we now begin to evaluate in the 
context of shadowing. 

In Figure [71 we illustrate the ability of this trans- 
port scheme to produce homogeneously illuminated radia- 
tion cones with an opening angle equal to the angular res- 
olution adopted for the scheme. In this example, a single 
source was placed in the center of a 2D unit square, and an- 
gular space has been divided into eight equal sized regions, 
with the source radiation only injected into four of them, al- 
ternating between an "empty" and a "full" cone. The green 
lines in the plot show some of the geometric boundaries of 
the angular discretisation as seen from the source. We see 
that the cone transport succeeds in producing a flat inten- 
sity proflle as a function of angle within every illuminated 
cone, while at the same time the leaking of radiation into 
cones that should remain dark as seen from the source is 
very small. We note that if we let the source inject radiation 
into two adjacent cones with equal luminosity, the radiation 



field shows no trace of the angular boundary between the 
cones, thanks to our use of the total intensity field in calcu- 
lating the local advection direction for the radiation of each 
partial field. 

An interesting question now arises how this transport 
scheme deals with obstacles and the problem of shadowing. 
We illustrate the salient points with a few tests in Figure [S] 
Here, we illuminate an obstacle (shown in white) by a single 
source in the left part of the simulated 2D space. We vary 
both the angular and the spatial mesh resolution in order 
to study how the shadowing performs in the cone transport 
scheme. In the left panel, we have used 50^ cell and eight 
angular regions. The fundamental cone size is shown by the 
green lines. In this setup, the opening angle of the obstacle 
as seen from the source is hence smaller than the angular res- 
olution of the RT, making the obstacle "unresolved" . In this 
case, the obstacle absorbs the correct amount of radiation 
expected for its size, but it will not form a correct shadow 
behind it. Instead, the "downstream region" behind the ob- 
stacle will get refilled with photons. As this can happen only 
by photons transported within the same geometric cone, a 
partial shadow is formed behind the object, with boundaries 
that are in principle parallel to the cone boundaries. In the 
middle panel, we repeat the test with the same spatial mesh 
resolution, but we have increased the number of cones to 32. 
In this way the angular size of the obstacle as seen from the 
source becomes larger than the angular resolution, allowing 
it to be resolved. As a result, a complete shadow is being 
formed, but this shadow is in general a bit smaller than the 
correct geometric shadow, with the difference being filled by 
a partial shadow, created in the cones that are only partially 
obscured by the obstacle. Finally, in the right hand panel, 
we have repeated the test on the left a second time, but now 
doubling the spatial resolution to 100^ cells while the angu- 
lar resolution was kept unchanged. The primary difference 
this makes is that the borders of the partial shadow that is 
formed are now more sharply defined compared to the case 
with lower spatial resolution, as expected. 

Finally, we examine how well our transport scheme can 
cope with different mesh geometries, which naturally arise 
in simulations with the AREPO code. In Figure [9] we show 
the radiation fields developing around a point source em- 
bedded in different mesh geometries: a Cartesian mesh, a 
hexagonal mesh (which is akin to the mesh geometry devel- 
oping in AREPO in regions of constant resolution), and an 
azimuthal/unstructured mesh. For all four cases, we com- 
pare the created radiation fields to the expected profile in 
2D, obtaining good agreement. This confirms the ability of 
our approach to work well with the unstructured Voronoi 
meshes produced by the AREPO code. 

3.2 Isothermal ionised sphere expansion 

We now turn to a test of our basic radiation advection 
scheme that involves both sources and sinks. To this end, 
we perform an ionised sphere expansion test in three dimen- 
sions, which is arguably the most fundamental and impor- 
tant test relevant for cosmic reionisation codes. 

The expansion of an ionisation front in a static, homo- 
geneous and isothermal gas is the only problem in radiation 
hydrodynamics that has a known analytical solution and is 
therefore indeed the most widely used test for RT codes. 
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Figure 8. Maps of the photon density field obtained with our "cone transport" sclieme for a point source at coordinates (0.1,0.5) and 
an obstacle (shown in white) centered at coordinates (0.72,0.5). The three panels differ in the mesh resolution used and the angular 
resolution employed for the radiation transfer. In all panels, the angular size of the cones employed in the angular discretisation are 
shown with green lines. In the left panel, where eight cones are used, the obstacle's opening angle as seen from the source is smaller 
than the fundamental cone, and therefore no complete shadow is formed. In the middle panel, 32 cones are used instead, such that the 
obstacle's opening angle is now larger than the angular resolution, allowing a full shadow to be formed. Finally, in the right hand panel, 
the spatial mesh resolution has been doubled in each dimension compared to the panel on the left, while the number of angular resolution 
elements has been kept at eight. Again there is no complete shadow formed, as expected, but the boundary of the shadow region behind 
the obstacle is now more sharply defined. 





Figure 9. Maps of the photon distribution for 2D simulations with three different cell shapes: Cartesian square, hexagonal and az- 
imuthal/unstructured. The line plot on the right shows the photon intensity profile, overplotted with the expected law. The vertical 
line indicates the average cell size. Results from all mesh shapes agree well with the analytical prediction (dot-dashed line). 



We adopt a monochromatic source that steadily emits A^.^ 
photons with energy hi/ = 13.6 eV per second into an ini- 
tially neutral medium with constant gas density riu- Then 
the Stromgren radius at which the ionised sphere around the 
source reaches its maximum radius is defined as 



rs,o 



1/3 



(32) 



where ob is the recombination coefficient. If we approximate 
the I-front is infinitely thin, i.e. features a discontinuity in 
the ionisation fraction, then the temporal expansion of the 
Stromgren radius can be solved analytically in closed form, 
with the I-front radius ri,o given by 



n,o = rs,o[l - exp{-t/tr. 



where 



il/3 



(33) 



(34) 



is the recombination time and qb is the recombination co- 
efficient. 

More accurately, the neutral and ionised fraction as a 
function of radius of the Stromgren sphere can be calcu- 
lated analytically (e.g. lOsterbrock fc Ferlandll2006l ) from the 
equation 

'^^^^^ J du N^{i/) e"^"^''^ a„ = fi^ii^r) nn as , (35) 

where nm is the neutral fraction, hmi is the ionised fraction 
and 



r^(r) = riHO-^ / dr' hm{r') 



(36) 



Moreover, we can analytically solve for the radial profile of 
the photon density n.y{r), yielding 



n (r)-i^ 



exp 



«;(r') dr' 



(37) 



From this we can also obtain the profile of the ionised 
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Figure 10. Profiles of the neutral and ionised fraction at the end 
of the ionised bubble expansion test, when the Stromgren radius 
rg has been reached. The black line is the analytic solution based 
on equation lj35|l , while the coloured lines are the numerical results 
for mesh resolutions of 20'^, 40'^, 80'^, and 160'^ cells, as labelled. 

fraction nHii(r) as a function of time. We note that the 
Stromgren radius obtained by direct integration of equation 
l|35p differs from the approximate expression p2l) because it 
does not approximate the ionised region as a top-hat sphere 
with constant ionised fraction. 

For definiteness, we follow in our tests the expansion 
of the ionised region around a source that emits A^.^ = 
5 X 10** s~^ photons. The surrounding hydrogen number 
density is set to wh ~ 10"'^ cm~^ at a temperature of 
T = 10'* K. At this adopted temperature, the case B recom- 
bination coefficient is Ob = 2.59x10"^"^ cm^s~^. Given these 
parameters, the recombination time is tree ~ 125.127 Myr 
and the expected Stromgren radius is rs,o = 5.38 kpc. 

In Figure IIOI we show the profiles of ionised and neu- 
tral fraction at the end of the ionised sphere expansion, when 
the Stromgren radius has been reached. We present results 
for simulations with four different spatial resolutions, using 
grids with 20^, 40^, 80^ and 160^ cells, respectively, using 
our point-source advection scheme. The results for all reso- 
lutions agree well with the analytical solution. The largest 
errors occur close to the central point source, but with bet- 
ter spatial resolution they become progressively smaller. In 
Figure [TTl we show the time evolution of the ionising front, 
for the same simulations. The position of the front is deter- 
mined as the distance from the source at which the ionised 
fraction equals 0.5. The agreement with the analytical solu- 
tion is generally good and improves with better resolution. 
However, in the beginning, the ionisation front moves no- 
ticeably slower than expected, which is due to our use of the 
reduced speed of light approximation with c' = c/1000. At 
later times, this initial error becomes unimportant, however, 
and the numerical solution matches the analytic expectation 
well. Making the start-up error vanishingly small would be 
possible, if desired, but requires using c' — c. 
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Figure 11. Radius of the ionised region as a function of time, 
in units of the recombination time tree- The solid black line is 
the analytic solution from equation I I35II . while the dashed black 
line gives the simple approximation of equation lj33p . The coloured 
lines give our numerical results for different mesh resolutions equal 
to 20^, 40^, 80^, and 160^ cells. We see that the ionising front 
is slower than the theoretical prediction in the beginning of the 
expansion, as a result of the artificially reduced speed of light 
adopted here. At late times, the numerical result agrees however 
very well with the analytical solution. 



0. 




Figure 12. Map of the neutral fraction in a slice through the 
center of our Stromgren sphere test (based on our point-source 
advection scheme), for our highest resolution simulation with 160^ 
mesh cells. The white line shows the contour at a neutral fraction 
of 0.5. 

In Figure [121 we show a map of the neutral fraction in a 
slice through the source plane for the resolution 160"^ . We no- 
tice that the isophotal shapes exhibit small departures from 
a perfectly spherical shape, which originate in spatial dis- 
cretisation errors close to the source. In fact, these deviations 
depend on the geometry of the source cell itself. For a hexag- 
onal mesh structure as it occurs for a regularised Voronoi 
mesh in 2D dimensions, the errors are noticeably smaller 
than for the Cartesian mesh employed here. Higher spatial 
resolution alone will normally not be able to decrease the de- 
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Figure 13. Neutral fraction in a slice through the center of 
two nearby sources of equal luminosity that ionise neutral gas. 
The three panels from top to bottom show different evolutionary 
stages. The top panel shows a stage before the ionised spheres 
overlap (t = 25Myr). Here they have exactly the same shape 
and do not influence each other yet. In the middle panel, the two 
have begun to overlap [t = 100 Myr), while in the bottom panel 
the final state is shown, where the ionised region becomes time 
invariant {t = 500 Myr). 




Figure 14. Neutral fraction along a line passing through the 
centres of two nearby sources that ionise the background gas. 
The green line shows the numerical result, whereas the black lines 
are a simple composite model for the expected structure of the 
solution based on superposing the analytic solution for each of the 
sources (gray dashed line for the left source and gray solid line 
for the right source) . This superposition of the individual sources 
describes the numerical solution reasonably accurately, but we 
note that it is not the correct solution; the latter can only be 
obtained numerically for this problem. 



0. 




Figure 15. Map of the neutral fraction in a slice through the 
center of a Stromgren sphere. Unlike in our previous tests, an 
absorbing obstacle in the form of an optically thick disk was in- 
cluded as well (shown as a black line) . We find that a nice shadow 
is produced behind the disk, with the inobscured directions de- 
veloping as in the Stromgren sphere without an obstacle. 

viations to arbitrarily small levels, but spreading the point 
source over multiple cells (effectively resolving the source 
geometry) can make the isophots perfectly round if desired. 
We note that our cone transport scheme also does a good 
job in producing round isophots, even when a single cell is 
used as source. 

As a simple variant of the isolated source case, we have 
also considered the evolution of the ionised regions around 
two sources that are 4kpc apart, using our multiple point- 
source scheme. The density of the gas and the luminosity of 
each source are the same as in the previous test. In Fig- 
ure 1131 we show maps of the neutral fraction in a slice 
through the source at three different times: t — 25 Myr 
(top), t = 100 Myr (middle), and t = 500 Myr (bottom). 
An important point of this test is that the proximity of the 
sources does not affect the shape of the ionised regions at 
all until they begin to overlap. This is very different in the 
OTVET scheme, for example, where the early expansion is 
distorted because the Eddington tensor estimates already 
"feel" nearby sources even though they may still be com- 
pletely hidden in their own ionisation bubbles. In Figure [T4l 
we show the neutral fraction along a line passing through 
both sources at the final time. A simple model for the ex- 
pected neutral fraction based on the superposition of the 
analytic single source solution is shown in black, while the 
numerical solution is shown in green. While the superpo- 
sition model does a reasonably good job in describing the 
numerical solution, we note that the latter is showing im- 
portant differences, for example for the radiation intensity 
between the sources. Our method allows an accurate calcu- 
lation of this quantity, and similarly for more complicated 
setups. 

In Figure [15] we show a further map of the neutral frac- 
tion in a slice through the source plane in a simple single- 
source Stromgren test. However, in this test we included 
an obstacle in the form of an optically-thick three dimen- 
sional plate, located 2kpc from the source (shown in black 
in the figure). The setup is meant to test shadowing in 3D 
for a problem with non-trivial source function, and is de- 
signed to match the p arameters of an equivalent test in 
iPawUk fc Schavd (|2008l ). We can see that our obstacle pro- 
duces a clear shadow that remains fully neutral, as expected. 
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Figure 16. Stromgren sphere test for our 'cone transport' scheme 
where the angular space is decomposed into cones of equal solid 
angle. The panel on top shows the profiles of ionised and neutral 
fractions (green lines) versus distance from the source, in units 
of the Stromgren radius rs q at the end of the expansion at t = 
500 Myr. The analytical solution is given by the solid black line. 
The bottom panel shows the ionisation front radius as a function 
of time, relative to the Stromgren radius rs g as a function of time 
in units of the recombination time tree ■ The numerical solution is 
shown in green and the analytical one in black. 



Comparing our result to those o f lPawlik fc Schavd (l2008l . see 
their Fig. 10), we find good qualitative agreement but much 
reduced numerical noise. 

Finally, we check whether using the cone transport 
scheme described in Section 12.41 is equally well capable of 
accurately solving the Stromgren sphere problem. To this 
end we have repeated our standard setup for the ionised 
sphere expansion of a single source, but this time employing 
direct discretisation of angular space using 12 cones for the 
full 47r solid angle, and a spatial mesh resolution of 40"^ . In 
the top and bottom panels of Figure [16] we show the pro- 
files of ionised and neutral fraction at the end of the ionised 
sphere expansion, and the temporal evolution of the ionising 
front, respectively. The numerical results agree well with the 
analytical solutions, with an overall accuracy that is com- 
parable to that of our point source treatment. 

3.3 Ionising front trapping in a dense clump 

In our next test, we study the behavior of the code in a 
more ch allenging se t ting t aken from the RT code comparison 
study of llliev et al.l (|2006l ). A plane-parallel front of ionising 
radiation is incident on a dense, cold clump. The I-front 
penetrates the clump, ionises it and heats it up. Eventually, 
the I-front gets trapped half-way through the clump, and as 
the it is stopped inside the obstacle, a shadow is produced 
behind the clump. 

Our set-up of this test problem is as follows. 




Figure 17. Maps of neutral fraction (top row) and temperature 
(bottom row) , in a simulation of the interaction of a plane-parallel 
ionisation front with a dense clump. The two columns show our 
simulation results at two different times, t = 1 Myr (left) and 
t = 15 Myr (right). We note that already at the earlier time the 
background gas has been fully ionised. The I-front gets however 
stuck in the clump, producing a shadow behind it. 
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Figure 18. Neutral fraction and temperature as a function of 
distance from the center of the dense clump, at three different 
times: t = I Myr (left column), 3 Myr (middle) and 15 Myr (right). 
The shaded area shows the geometric extension of the clump. Re- 
sults obtained with the code CRASH in the RT code comparison 
project are also included and shown as dashed lines. 



We simulate a plane-parallel I-front with fiux F-y — 
lO'' photons s~^ cm"'^ that is incident on a dense clump, 
located 5 kpc away from the edge of the simulation do- 
main. The ambient background gas has density nn = 2 x 
10"" cm"^ and temperature T = 8000 K. The radius of the 
clump is Tciump = 0.8 kpc, with a density of n^"™'' = 200 nn, 
and a te mperature rdump = 40 K. We note that in this test, 
following llliev et aL I l|2006l ), the gas is not allowed to move 
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3.4 lonisation of a static cosmological density field 

In our most demanding test of pure RT we follow hydro- 
gen lonisation in a realistic cosmological density field, which 
is taken to be static for simplicity. Again, in order to be 
able to compare our results with those of th e cosmolog- 
ical RT comparison project jlliev et alj l2006h we use the 
same density field, the same cosmological box parameters, 
and assign sources in the same way. The test is based on 
a cosmological density field in a periodic box with size 
0.5 ft^^comoving Mpc that resulted from the evolution of a 
standard ACDM model with the cosmological parameters 
no = 0.27, r^b = 0.043, and h = 0.7. The gas density field 
at redshift 2: = 9 is considered for further analysis. 

The source distribution is determined by finding halos 
within the simulation box with a Friend-Of-Friends (FOF) 
algorithm and then assigning sources to the 16 most massive 
groups. The photon luminosity of these sources is taken to 
be 



Figure 19. Time evolution of the temperature, neutral fraction, 
and ionised front (solid lines) in the dense clump that is ionised 
by an impinging I-front. We also include results obtained with 
the code CRASH in the RT code comparison project, which are 
shown as dashed lines for comparison. 



due to pressure or gravitational forces, hence only radia- 
tive transfer is tested. The system is evolved for a period of 
15 Myr with a resolution of 40^ cells. In Figure 1171 we show 
the neutral fraction and the temperature of the system in 
slices through the centre of the clump at times t — 1 Myr, 
and 15 Myr. The I-front approaches from the left, moving 
to the right. At time t = 1 Myr, already the whole box has 
been swept up by the ionising photons, with the clump pro- 
ducing a clear shadow in the downstream direction on the 
right hand side of the clump. As time advances further, the 
clump becomes more ionised and continues to heat up, but 
the shadow is preserved throughout the simulated time span 
without being filled in by diffusion. 

In Figure [18] we show the temperature and neutral frac- 
tion as a function of distance from the geometric clump 
center. The results are compared to those obtained in the 
comparison study of llliev et al.l l|2006l ) for the Monte-Carlos 
transfer code CRASH. The position and shape of the ionis- 
ing front agree well with the results from the CRASH code, 
both at times t = 1 Myr and 15 Myr. The temperature pro- 
file shows, however, some differences. This discrepancy can 
be traced back to inaccuracies in CRASH, where lower en- 
ergy photons penetrate into the gas ahead of the ionising 
front and heat it there. 

Finally, in Figure [191 we show the time evolution of the 
temperature, ionised fraction and position of the I-front in 
the clump, compared to the results obtained with CRASH. 
The clump is 60% ionised at the end, its average temperature 
increases to several lO^K and the I-front becomes trapped 
around the geometric center of the clump, which is all in 
good agreement with the CRASH results. We hence conclude 
that our RT scheme yields results of good accuracy for this 
test, which are comparable in accuracy to those obtained 
with expensive yet accurate Monte-Carlo treatments. 



norupts 



(38) 



where M is the total halo mass, ts = 3 Myr is the assumed 
lifetime of each source, rrip is the proton mass, and f-y = 
250 is the number of emitted photons per atom during the 
lifetime of the source. For simplicity we also set the initial 
temperature of the gas to 100 K throughout the whole box. 

In Figure [20I we show the neutral fraction and the tem- 
perature in slices through the center of the simulated vol- 
ume, at the final evolution time of t = 0.4 Myr. We have 
calculated the radiation transfer in three different ways, cor- 
responding to the three variants of the radiation advection 
approach proposed in this paper. In the left panel, we show 
the results if all sources are treated as linearly independent 
point sources. The middle panel shows the result when only 
the Abr — 4 brightest sources seen from a given point are 
treated as point sources, while the remaining luminosity is 
fed to a background radiation field which is treated with 
radiative diffusion. Finally, the right hand panel gives the 
result when angular discretisation with 12 HEALPIX cones 
for the full solid angle is applied. 

Visually, based on Fig. 1201 all three results agree very 
well with each other. However, there are some small differ- 
ences in the structure of the ionised regions and in the shape 
of the I-fronts. For a better quantitative comparison we show 
in the top panel of Figure [21] the volume filling function of 
the neutral fractions for all three simulation methods, where 
a comparison with CRASH res ults from the RT code com- 
parison study (|lhev et al.ll200^ ) is also included. Our results 
agree very well with the CRASH data, and we note that there 
are also no substantial differences between the three different 
approaches for dealing with multiple sources in our radia- 
tion advection scheme. The same conclusion is also reached 
from a comparison of the volume distribution function of 
the temperature, which is shown in the bottom panel of Fig- 
ure 1211 We note that these results are considerably better 
than those we obtained for t he OTVET scheme implem ented 
in the SPH code GADGET (jPetkova fc Springel|[2009l ). 
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Figure 20. Maps of the ionised fraction (top row) and temperature (bottom row) in a slice through the middle of the simulation volume 
at time t = 0.4 Myr of our cosmological density field ionisation test. In the results shown in the left column, all sources have been 
treated independently as point sources. In the middle column, only the locally four brightest sources have been considered independently, 
while the remaining luminosity has been treated with radiation diffusion. Finally, the results in the right column are based on our cone 
transport algorithm with a division of the full solid angle into 12 cones of equal size, corresponding to the coarsest HEALPIX resolution. 




Figure 21. Distribution functions of the neutral fraction (top panel) and the temperature (bottom panel) in the simulated cosmological 
volume at time t = 0.4 Myr, for three different variants of our radiation transfer scheme, as labelled: (1) all sources are followed in 
a linearly independent fashion, (2) only the four locally brightest sources are followed as point sources with the rest treated through 
radiative diffusion, and (3) a cone transport approach based on a division of the u nit sphere into 12 cones. For comparison, we also 
include results obtained with the code CRASH in the RT code comparison project of llliev et al. I ||2006| '). 
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Figure 22. Profiles of ionised fraction, hydrogen number density, pressure, temperature, and mach number at different times for a 
hydrodynamically coupled Stromgren sphere test. The distance from the source is normalised by the box size L^ox = 15kpc. The three 
lines in each plot correspond to the times t = lOMyr (solid), t = 200 Myr (dashed), and t = 500 Myr (dotted). 



3.5 Ionised sphere expansion in a dynamic 
density field 

As our final test, we again follow the expansion of an ionised 
sphere in an initially homogeneous and isothermal medium, 
similar to Section [3.21 but this time we allow the gas to be 
heated up by the photons and to expand due to the raised 
pressure. This is hence a radiation hydrodynamics test where 
both RT and hydrodynamics are f ollowed. We desig n our 
test similar to the set-up studied in llliev et al. I l|2009l ). The 
source is at the center of the simulation domain and emits 
at a luminosity of Nj = 5 x 10'** s~^ . The surrounding 
hydrogen number density is nu = 10~^ cm~^ at an initial 
temperature of T = 10^ K. The simulated box is 30kpc on 
a side, and is resolved with 160^ cells. We evolve the system 
for 500 Myr. 

Th ere are two c ritical gas velocities defined for such a 
set up (|Spitzerlll978l ): the R-critical velocity v-r = and 
the D-critical velocity vo ~ ~ ^ (4)^ — (cj)^, where 
and Cg are the isothermal sound speeds ahead and behind 
the I-front, respectively. If we assume the ionised gas has 
temperature 10'* K and the neutral gas 10^ K, than we obtain 

~ 25.70 kms~^ and ud ~ O.OSkms"'^. The I-front is 
called D-type when its speed is smaller than the D-critical 
speed vi ^ ud • In this case it is subsonic with respect to the 
neutral gas, which expands as the I-front passes through it. 
When vi ^ Vr, the I-front is called R-type. It is supersonic 
with respect to the neutral gas ahead, and the gas does not 



expand as the I-front passes trough it. When wd < wi < vr, 
there is a hydrodynamic shock wave in front of the I-front. 
The position of the I-front in this stage is given by l|Spitzeij 

[Hzi) 

where Cs is the sound speed of the ionised gas and rs,o is the 
Stromgren radius given by equation ((32]). 

We evolve our test setup for 500 Myr and analyse the 
results at three different times equal to f = 10, 200 and 
500 Myr. Figure [22] shows the time evolution of the profiles 
of the ionised fraction, hydrogen number density, pressure, 
temperature, and mach number profiles. At time t = 10 Myr, 
the gas expands at subsonic speed. The pressure inside the 
ionised bubble is very high as the density is still relatively 
close to 10~^ cm~'^ and the temperature is several 10* K. At 
later times, t — 200 Myr, there is a shock developing ahead 
of the ionising front. The gas in this pseudo shock region is 
compressed, leading to densities higher than 10~^ cm~^ and 
an increased pressure. At time t = 500 Myr, the dynamic 
situation of the gas is similar, as there is still a shock ahead 
of the ionising front, but the pressure in the ionised bubble 
has dropped significantly due to the lowered density of the 
gas. 

In Figure [23] we show the evolution of the radius of the 
I-front. In the first 40 Myr, the ionising front moves with a 
speed larger than the R-critical velocity: vi > wr. Its evolu- 
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Figure 23. Evolution of the position of the ionising front in a 
hydrodynamically coupled Stromgren sphere test. The distance 
is expressed in terms of the Stromgren radius rs q for the case of 
a static density field. The dotted line shows the analytic solution 
for the time evolution in the static density case, while the dashed 
line gives the solution for the dynamic density case. The latter 
is well reproduced by our numerical AREPO calculation. In the 
bottom panel, we show the speed of the ionisation front. In the 
first 40 Myr of the expansion, the front moves with a speed higher 
than the R-critical velocity (indicated by a dotted line). 

tion corresponds to that of an I-front in a static density field, 
and the position of the front follows the analytical prediction 
from equation (|33p . The gas does not expand significantly 
in this stage. As the speed of the I-front drops below the 
R-critical velocity, a shock develops ahead of it and the gas 
gets compressed in these regions. Here the position of the 
front evolves according to equation H39[l . 

In general, our results for this t est agree w e ll wit h the 
other codes that have been tested bv llliev et al.l (|2009l '). We 
find the best agreement with the ENZO-RT results from the 
RT comparison study, which is probably due to the specific 
monochromatic nature of our code. 



4 DISCUSSION AND CONCLUSIONS 

In this study, we have proposed a novel implementation of 
radiative transfer and implemented it in the moving-mesh 
code AREPO. The method differs substantially from com- 
monly employed ray-tracing or moment-based schemes in 
that it directly evolves a discretised version of the Boltz- 
mann equation describing the photon distribution function. 
This is done in terms of an advection treatment, where the 
photon transport is carried out with a second-order accu- 
rate upwind scheme, based on methods that are commonly 
employed in hydrodynamic mesh codes. We have introduced 
three different approaches to deal with multiple sources, ei- 
ther by splitting up the radiation field into a linear sum 
of the partial fields created by all sources, by using a hy- 
brid approach consisting of an exact treatment of the lo- 
cally brightest sources combined with radiative diffusion, or 
by employing a direct discretisation of angular space into a 



finite set of cones. The latter approach is the most general. 
At a given angular resolution, it can easily deal with an arbi- 
trary number of sources as well as with radiation scattering. 
Also, if the number of angular cones is enlarged, its angu- 
lar accuracy becomes progressively better, allowing a simple 
way to test for convergence with angular resolution. 

The radiation transport in our method is manifestly 
photon conserving. Combined with a photon-conserving 
treatment of the source terms, this yields a very robust de- 
scription of the reionisation problem, ensuring that ionisa- 
tion fronts propagate at the correct speed. If needed, our 
code can employ a reduced speed-of-light approximation 
that avoids overly small timesteps while not altering the 
growth of ionised regions in any significant way. 

We have presented tests of our new scheme in a vari- 
ety of cases. Using different photon transport tests in 2D, 
we have shown that our method manages to accurately cap- 
ture shadows, and to produce the correct radiation fields 
independent of the mesh geometry. To test the coupling 
of gas physics with the photon transport we have carried 
out isothermal ionised sphere expansion tests and compared 
to the analytical solutions. The results agree reassuringly 
well with theoretical expectations, both for our linear sum- 
mation method and the cone transport approach. Further- 
more, we have shown that our method can treat multiple 
point sources in a highly accurate way, without the prob- 
lem of a detrimental mutual influence of the sources onto 
each othe r, which is encounter e d in certain mom ent-based 
schemes l|Gnedin fc Abell l200ll : iPetkova fc SprTn gcl 200^). 
We have also shown that our code performs well on the prob- 
lem of the ionisation of a static cosmological density field, 
where we benchmarked our results against those obtained for 
th e same setu p in the radiative transfer comparison study 
of llliev et af] |2006y Similarly, our results for I-front trap- 
ping in a dense clump, and for a hydrodynamically coupled 
Stromgren sphere agree well with th ose of other radiative 
transfer codes included in llliev et al. I iiooe). 

Compared to other radiative transfer schemes, our new 
method based on the cone transport features several in- 
teresting advantages. Unlike long-characteristics or Monte 
Carlo schemes, it avoids any strong sensitivity of the com- 
putational cost on the number of sources, and it does not 
concentrate the computational effort in regions close to the 
sources, which greatly helps in parallelising the calculations. 
Also, the ability to easily treat time-dependent effects that 
are consistently coupled to the hydrodynamic evolution is a 
substantial asset. Compared to moment-based solvers, our 
method can cast sharp shadows, and it performs accurately 
both in the optically thin and the optically thick regimes. 

Our cone-transport scheme b ears some superficial re - 
sembles to the TRAPHIC scheme of lPawlik fc Schavel (|2008l) . 
However, unlike their approach, we do not rely on stochastic 
Monte Carlo techniques. Instead, we work with an explicit 
spatial reconstruction of the radiation field and a fixed set 
of angular cones. As a result, our radiation field is essen- 
tially free of stochastic noise, which is a significant advan- 
tage compared to Monte Carlo approaches and offers much 
better convergence rate. 

We hence think that our new method represents a 
promising treatment of radiative transfer in cosmological 
simulations. Its simple and general principles make it appli- 
cable in a wide variety of different problems in astrophysics. 
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In particular, it should not only be useful for studying reion- 
isation of the Universe, but also, for example, for studying 
star formation in molecular clouds, wher e ionising radiation 
creates pillars or neutral dense gas (e.g. iGritschneder et aU 
l2009t ). Also, our RT solver coupled to AREPO should al- 
low self-consistent and more accurate treatments of radiative 
feedback effects in hydrodynamic simulations of star forma- 
tion, something that we will study in forthcoming work. 
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